function [dets, suminv] = scgibbshelper(a,b,K)
% X_k = a*I_{k*k} + b*1_{k*k}
% e.g. [[a+b b b b];
%       [b a+b b b];
%       [b b a+b b];
%       [b b b a+b]];

% a, b:   matrices of same dimension
% dets:   dets(k) = det(X_k)
% suminv: suminv(k) = sum(sum(inv(X_k)))

if isvector(a) && isvector(b)
    a = a(:); b = b(:);
    dets = zeros(numel(a), K);
    suminv = zeros(size(dets));
    dets(:,1:2) = [a+b, (a+2*b).*a];
    suminv(:,1:2) = [1./(a+b), 2./(a+2*b)];

    apkm1 = a;
    for k = 3:K
        apkm1 = apkm1.*a;
        tmp = k.*b + a;
        dets(:,k) = apkm1.*tmp;
        suminv(:,k) = k./tmp;
    end
else
    dets = zeros([size(a), K]);
    suminv = zeros(size(dets));
    dspan = cell(1, ndims(a));
    dspan(:) = {':'};
    dets(dspan{:},1) = a+b;
    dets(dspan{:},2) = (a+2*b).*a;
    suminv(dspan{:},1) = 1./(a+b);
    suminv(dspan{:},2) = 2./(a+2*b);

    apkm1 = a;
    for k = 3:K
        apkm1 = apkm1.*a;
        tmp = k.*b + a;
        dets(dspan{:},k) = apkm1.*tmp;
        suminv(dspan{:},k) = k./tmp;
    end
end

end

% % Numerical Stability Test
% b = 5e8*rand(100000, 1);
% a = ones(size(b));
% [dets, suminv] = scgibbshelper(a, b, 5);
% newk = 3;
% inva = dets(:,newk-1)./dets(:,newk);
% invb = -inva+0.5*dets(:,2).*suminv(:,2).*dets(:,newk-2)./dets(:,newk);
% tSigma = repmat(shiftdim(invb',-1), [newk, newk, 1]);
% for kk = 0:newk-1
%     tSigma(kk*(newk+1)+1:newk*newk:end) = inva';
% end
% normerr = zeros(1,size(b,1));
% for i = 1:size(b,1)
%     M = a(i)*eye(newk)+b(i)*ones(newk);
%     normerr(i) = norm(inv(M)-tSigma(:,:,i), 'fro');
%     if det(tSigma(:,:,i)) < 0
%         fprintf('%d,', i);
%     end
% end
% fprintf('\n');